Trabajo final del Diplomado en AI & DL Universidad Nacional de Colombia

Clasificación de cambios de coberturas vegetales en imágenes de satélite Sentinel-2 en Bogotá

Fecha: Mayo de 2021

Autores: Samuel Mesa samuelmesa@linuxmail.org, Wilgen Correa

Zona de estudio: Tunjo

Cargar librerias

import datetime
import os, glob
import numpy as np
import seaborn as sns
import geopandas as gpd
import pandas as pd
import imageio
from pyspatialml import Raster
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import tensorflow as tf
from tensorflow.keras.layers import Dense, Dropout
from tensorflow.keras.models import Sequential
from sklearn.impute import SimpleImputer
from tensorflow.keras.wrappers.scikit_learn import KerasClassifier
from tensorflow.keras.callbacks import CSVLogger, RemoteMonitor
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import accuracy_score, precision_score, recall_score, confusion_matrix
import earthpy.plot as ep
import earthpy.spatial as es

Realizar el cargue de imágenes y las áreas de enetrenamiento

Bandas a seleccionar de la imagen de entrada de SENTINEL-2, con remuestreo a 10 m de resolución espacial

1. B2 (Blue)
2. B3 (Green)
3. B4 (Red)
4. B5 (Veg. red edge)
5. B6 (Veg. red edge)
6. B7 (Veg. red edge)
7. B8 (NIR)
8. B8A (Narrow NIR)
9. B11 (SWIR)
10. B12 (SWIR)

La estructura de directorios es el siguiente.

+-- [nombre_proyecto]

+-- shapes: shapes de aoi y capas de superposición

+-- models: Almacena el modelo

+-- sources: imágenes satelitales de entrada

+-- results: resultados de las clasificación

+-- logs: logs del entrenamiento

+-- figures: imágenes y figuras

path_project = "./tunjo/"
path_sources = os.path.join(path_project, "sources")
path_shapes = os.path.join(path_project, "shapes")
path_results = os.path.join(path_project, "results")
path_logs = os.path.join(path_project, "logs")
path_figures = os.path.join(path_project, "figures")
path_models = os.path.join(path_project, "models")
list_paths = [path_results, path_figures, path_models, path_logs]

# Crear directorios temporales

for path in list_paths:  
    try:
        os.mkdir(path)
    except FileExistsError:
        print("Directorio ya existe: ", os.path.basename(path))

# Defina la imagen de entrada 
img_train = '20210309T152639_20210309T152638_T18NWL.tif'
img_file = os.path.join(path_sources, img_train) 
aoi_file = os.path.join(path_shapes, 'aoi.geojson') # Formato geográfico
manzana_file = os.path.join(path_shapes, 'manzana.geojson') # Formato geográfico
img_name = img_train.split('_')[0] 

# Cargar la imagen como un objeto Raster Dataframe y el shape de aoi como Geodaataframe
stack = Raster(img_file)
training = gpd.read_file(aoi_file)
manzana = gpd.read_file(manzana_file)

bandsio = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12']
bandsio_names = ['B2 (Blue)', 'B3 (Green)', 'B4 (Red)', 'B5 (Veg. red edge)', 
               'B6 (Veg. red edge)', 'B7 (Veg. red edge)', 'B8 (NIR)', 
               'B8A (Narrow NIR)', 'B11 (SWIR)', 'B12 (SWIR)']

# Mapea los nombres de bandas de entrada
[stack.rename({name: bandsio[idx]}, in_place=True) for idx, name in enumerate(stack.names)]
extent = es.rio.plot.plotting_extent(stack)

Directorio ya existe:  results
Directorio ya existe:  figures
Directorio ya existe:  models
Directorio ya existe:  logs

Se define la selección de bandas de entrenamineto y salida para todo el proceso

bandsout = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12']
bandsout_names = ['B2 (Blue)', 'B3 (Green)', 'B4 (Red)', 'B5 (Veg. red edge)', 
               'B6 (Veg. red edge)', 'B7 (Veg. red edge)', 'B8 (NIR)', 
               'B8A (Narrow NIR)', 'B11 (SWIR)', 'B12 (SWIR)']

Anáisis exploratorio de datos

Despliegue de imágenes y análisis exploratorio

# Despliegue los canales  imágenes
ep.plot_bands(stack[bandsout].read(), title=bandsout_names, figsize=(10, 8))

array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7ffa4462d730>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c62a3d90>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c621dcd0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c61bd130>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c61e1550>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c618aa30>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c61a8f10>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c61b4fa0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c6111700>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c60bbb20>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c60e6eb0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9c609d3d0>]],
      dtype=object)

# Histogramas por bandas

colors_list = ["Blue", "Green", "Red", "Salmon", "Tomato", "Coral", "Orangered",
               "Chocolate","Darkorange","Maroon"]

ep.hist(stack[bandsout].read(), colors=colors_list, title=bandsout_names, 
        ylabel='Reflectancia', bins=50, cols=2)

(<Figure size 864x864 with 10 Axes>,
 array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bd6e3640>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bd70b0d0>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc6c57c0>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc671b20>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc61fe80>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc5d82e0>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc5f6670>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc6026d0>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc55bca0>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x7ff9bc5870a0>]],
       dtype=object))

# Despliegue composición de color

rgb432 = stack[['B4','B3', 'B2']].read()
rgb843 = stack[['B8','B4', 'B3']].read()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

ep.plot_rgb(rgb432, ax=ax1, title='Despliegue verdadero color (RGB432)', stretch=True)
ep.plot_rgb(rgb843, ax=ax2, title='Despliegue falso color (RGB843)', stretch=True)
plt.tight_layout()
plt.show()

Despliegue de las áreas de entrenamiento

fig, ax = plt.subplots(figsize=(9, 9))
ep.plot_rgb(rgb843, ax=ax, stretch=True, extent=extent)
manzana.boundary.plot(ax=ax, color='grey', alpha=0.5,)
training.plot(column="label", cmap='RdYlGn', ax=ax, legend=True,  alpha=0.65, categorical=True)
ax.axis('on')
plt.title('Áreas de entrenamiento')
plt.show()

# Convertir los pixeles a un dataframe de pandas

df = stack[bandsout].to_pandas()
df.head()

	x	y	B2	B3	B4	B5	B6	B7	B8	B8A	B11	B12
0	594300.000000	506180.0	306.0	552.0	521.0	764.0	1610.0	2268.0	2382.0	2557.0	1628.0	1033.0
1	594310.068493	506180.0	248.0	354.0	308.0	764.0	1610.0	2268.0	1790.0	2557.0	1628.0	1033.0
2	594320.136986	506180.0	290.0	575.0	449.0	973.0	2427.0	2693.0	2808.0	2838.0	1856.0	1064.0
3	594330.205479	506180.0	334.0	666.0	514.0	973.0	2427.0	2693.0	3164.0	2838.0	1856.0	1064.0
4	594340.273973	506180.0	363.0	743.0	607.0	1132.0	2492.0	2963.0	2900.0	3366.0	2057.0	1234.0

print('Las columnas son: ', df.columns)
print('El tamaño´ del dataframe: ', df.shape)
print(df[bandsout].describe().T)

Las columnas son:  Index(['x', 'y', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11',
       'B12'],
      dtype='object')
El tamaño´ del dataframe:  (18816, 12)
       count         mean         std    min     25%     50%      75%     max
B2   18816.0   836.277583  530.597252   64.0   454.0   769.0  1080.00  5908.0
B3   18816.0  1066.075999  532.072826  133.0   722.0  1017.0  1284.00  6556.0
B4   18816.0  1138.359800  609.381740  114.0   714.0  1116.0  1456.00  6544.0
B5   18816.0  1406.380102  503.829649  308.0  1094.0  1391.0  1648.00  5772.0
B6   18816.0  1920.350872  486.753114  319.0  1608.0  1870.5  2184.25  5491.0
B7   18816.0  2132.852572  568.216038  326.0  1734.0  2036.0  2496.25  5425.0
B8   18816.0  2153.274979  655.596387  273.0  1678.0  2034.0  2602.00  6136.0
B8A  18816.0  2261.102891  617.403156  307.0  1804.0  2126.0  2699.00  5409.0
B11  18816.0  2173.782100  563.845698  392.0  1848.0  2155.5  2457.00  5675.0
B12  18816.0  1824.713010  675.126151  277.0  1379.0  1856.0  2210.00  6132.0

# Correlacción de las bandas

correlation_data = df.iloc[:,2:].corr()
correlation_data.style.background_gradient(cmap='coolwarm', axis=None)  

	B2	B3	B4	B5	B6	B7	B8	B8A	B11	B12
B2	1.000000	0.972844	0.952348	0.854090	0.278757	0.065329	0.073105	-0.026923	0.611698	0.808887
B3	0.972844	1.000000	0.971263	0.895550	0.395855	0.185262	0.199752	0.092389	0.688329	0.822943
B4	0.952348	0.971263	1.000000	0.910007	0.340398	0.118602	0.120648	0.022096	0.695865	0.843949
B5	0.854090	0.895550	0.910007	1.000000	0.557775	0.339255	0.275435	0.245003	0.833778	0.883065
B6	0.278757	0.395855	0.340398	0.557775	1.000000	0.948065	0.872684	0.913443	0.687255	0.358862
B7	0.065329	0.185262	0.118602	0.339255	0.948065	1.000000	0.915475	0.976660	0.538896	0.155824
B8	0.073105	0.199752	0.120648	0.275435	0.872684	0.915475	1.000000	0.917705	0.486189	0.119533
B8A	-0.026923	0.092389	0.022096	0.245003	0.913443	0.976660	0.917705	1.000000	0.502150	0.097535
B11	0.611698	0.688329	0.695865	0.833778	0.687255	0.538896	0.486189	0.502150	1.000000	0.868317
B12	0.808887	0.822943	0.843949	0.883065	0.358862	0.155824	0.119533	0.097535	0.868317	1.000000

# Extraer información de las bandas con la áreas de entrenamiento

df_shape = stack[bandsout].extract_vector(training)

df_shape = df_shape.merge(
    right=training.loc[:, ["label", "id"]],
    left_on="geometry_idx",
    right_on="index",
    right_index=True
)

df_shape = df_shape.dropna()

# Descripcion de los datos para las areas de entrenamineto

print('Total canales espectrales: ',  df_shape.columns)
print('Tamaño de entrenamiento: ', df_shape.shape)
print(df[bandsout].describe().T)

Total canales espectrales:  Index(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12',
       'geometry', 'label', 'id'],
      dtype='object')
Tamaño de entrenamiento:  (3370, 13)
       count         mean         std    min     25%     50%      75%     max
B2   18816.0   836.277583  530.597252   64.0   454.0   769.0  1080.00  5908.0
B3   18816.0  1066.075999  532.072826  133.0   722.0  1017.0  1284.00  6556.0
B4   18816.0  1138.359800  609.381740  114.0   714.0  1116.0  1456.00  6544.0
B5   18816.0  1406.380102  503.829649  308.0  1094.0  1391.0  1648.00  5772.0
B6   18816.0  1920.350872  486.753114  319.0  1608.0  1870.5  2184.25  5491.0
B7   18816.0  2132.852572  568.216038  326.0  1734.0  2036.0  2496.25  5425.0
B8   18816.0  2153.274979  655.596387  273.0  1678.0  2034.0  2602.00  6136.0
B8A  18816.0  2261.102891  617.403156  307.0  1804.0  2126.0  2699.00  5409.0
B11  18816.0  2173.782100  563.845698  392.0  1848.0  2155.5  2457.00  5675.0
B12  18816.0  1824.713010  675.126151  277.0  1379.0  1856.0  2210.00  6132.0

fig, axes = plt.subplots(1,2, figsize=(10,4), sharey=True, sharex=True)
plt.suptitle("Áreas de entrenamiento")
axes[0].set_title('No vegetación')
axes[1].set_title('Vegetación')
axes[0].set_xlabel('Bands')
axes[0].set_ylabel('Reflectance')
axes[1].set_xlabel('Bands')
sns.pointplot(data=df_shape[df_shape['id'] == 0][bandsout], ax=axes[0], 
              scale=0.5, estimator=np.mean)
sns.boxplot(data=df_shape[df_shape['id'] == 0][bandsout], ax=axes[0])
sns.pointplot(data=df_shape[df_shape['id'] == 1][bandsout], ax=axes[1], 
              scale=0.5, estimator=np.mean)
sns.boxplot(data=df_shape[df_shape['id'] == 1][bandsout], ax=axes[1])
axes[0].grid()
axes[1].grid()
plt.tight_layout()
fig.show()

fig, ax = plt.subplots(figsize=(10,5), sharey=True, sharex=True)
ax.set_title('Firmas espectrales de las áreas de entrenamiento')
ax.set_xlabel('Bands')
ax.set_ylabel('Reflectance')
sns.pointplot(data=df_shape[df_shape['id'] == 0][bandsout], ax=ax, 
              scale=1, estimator=np.mean, color='red')
sns.pointplot(data=df_shape[df_shape['id'] == 1][bandsout], ax=ax, 
              scale=1, estimator=np.mean, color='green')
ax.grid()
fig.show()

Entrenamineto de la red neuronal

Dividir las areas de entrenamineto y validación

#bands = ['B2', 'B3', 'B4', 'B7']

X = df_shape[bandsout].values
y = df_shape["id"].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

# Logs y callbacks para el entrenamiento

timename = datetime.datetime.now().strftime("%Y%m%d-%H%M%S")

log_csv = os.path.join(path_logs, f"log_{img_name}-{timename}.csv")
model_h5 = os.path.join(path_models, f"model_{img_name}-{timename}.h5")
imgkeras_classification = os.path.join(path_results, f"classkeras_{img_name}-{timename}.tif")
imgkeras_probability = os.path.join(path_results, f"probkeras_{img_name}-{timename}.tif")
imgkneig_classification = os.path.join(path_results, f"classkneig_{img_name}-{timename}.tif")
imgkneig_probability = os.path.join(path_results, f"probkneig_{img_name}-{timename}.tif")
fig_model = os.path.join(path_figures, f"train_{img_name}-{timename}.png")
movie_classkeras = os.path.join(path_figures, f"movie_classkeras_{timename}.gif")

call_save_model = tf.keras.callbacks.ModelCheckpoint(
    model_h5, monitor='val_loss', verbose=0, save_best_only=False,
    save_weights_only=False, mode='auto', save_freq='epoch', options=None
)

call_tensorboard = tf.keras.callbacks.TensorBoard(
    log_dir='./logs', histogram_freq=1, write_graph=True, write_images=True,
    update_freq='epoch', profile_batch=2, embeddings_freq=0,
    embeddings_metadata=True
)

call_csv = CSVLogger(log_csv, separator=",", append=False)

call_remote = RemoteMonitor(
    root='http://localhost:9000',
)

# Crear la funcion model de keras

def create_keras_model(layer1_units, layer2_units, dropout_rate, l2_regularization):
    model = Sequential([
        ## Adiciona capas una por una
        Dense(units=layer1_units, activation='relu', input_shape=(len(bandsout),)),
        # Adding dropout to prevent overfitting (regularización)
        Dropout(dropout_rate), # 10% out in each epoc
        Dense(units=layer2_units, activation='relu'),
        # Adding dropout to prevent overfitting (regularización)
        #model.add(Dropout(0.1))
        Dense(units=1, activation='sigmoid')
    ])
    
    model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
    print(model.summary())
    return model

# Realizar el entrenamiento

classifier_keras = KerasClassifier(
    build_fn=create_keras_model, 
    batch_size=32, 
    layer1_units=16,
    layer2_units=32,
    dropout_rate=0.025,
    l2_regularization=0,
    epochs=50, 
    shuffle=True, 
    validation_split=0.2,
    callbacks=[call_save_model, call_tensorboard, call_csv, call_remote],
    verbose=True
)

pipeline_keras = Pipeline([
    ('imputer', SimpleImputer()), 
    ('scaler', StandardScaler()), 
    ('model', classifier_keras)
])

pipeline_keras.fit(X=X_train, y=y_train)

Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense (Dense)                (None, 16)                176       
_________________________________________________________________
dropout (Dropout)            (None, 16)                0         
_________________________________________________________________
dense_1 (Dense)              (None, 32)                544       
_________________________________________________________________
dense_2 (Dense)              (None, 1)                 33        
=================================================================
Total params: 753
Trainable params: 753
Non-trainable params: 0
_________________________________________________________________
None
Epoch 1/50
68/68 [==============================] - 3s 26ms/step - loss: 0.6184 - accuracy: 0.6805 - val_loss: 0.2192 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 2/50
68/68 [==============================] - 0s 5ms/step - loss: 0.1466 - accuracy: 0.9985 - val_loss: 0.0409 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 3/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0310 - accuracy: 0.9998 - val_loss: 0.0153 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 4/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0140 - accuracy: 0.9984 - val_loss: 0.0086 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 5/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0098 - accuracy: 0.9992 - val_loss: 0.0054 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 6/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0060 - accuracy: 1.0000 - val_loss: 0.0040 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 7/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0039 - accuracy: 0.9997 - val_loss: 0.0032 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 8/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0035 - accuracy: 0.9999 - val_loss: 0.0025 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 9/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0021 - accuracy: 1.0000 - val_loss: 0.0022 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 10/50
68/68 [==============================] - 0s 6ms/step - loss: 0.0032 - accuracy: 0.9985 - val_loss: 0.0016 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 11/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0023 - accuracy: 0.9995 - val_loss: 0.0014 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 12/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0020 - accuracy: 0.9995 - val_loss: 0.0012 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 13/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0015 - accuracy: 0.9996 - val_loss: 0.0011 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 14/50
68/68 [==============================] - 0s 4ms/step - loss: 0.0016 - accuracy: 1.0000 - val_loss: 0.0010 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 15/50
68/68 [==============================] - 0s 4ms/step - loss: 8.3093e-04 - accuracy: 1.0000 - val_loss: 9.3394e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 16/50
68/68 [==============================] - 0s 4ms/step - loss: 9.8381e-04 - accuracy: 1.0000 - val_loss: 9.1729e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 17/50
68/68 [==============================] - 0s 5ms/step - loss: 8.8929e-04 - accuracy: 1.0000 - val_loss: 7.3446e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 18/50
68/68 [==============================] - 0s 4ms/step - loss: 9.0982e-04 - accuracy: 1.0000 - val_loss: 6.3256e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 19/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0012 - accuracy: 1.0000 - val_loss: 4.8234e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 20/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0014 - accuracy: 0.9990 - val_loss: 6.0310e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 21/50
68/68 [==============================] - 0s 5ms/step - loss: 0.0010 - accuracy: 1.0000 - val_loss: 5.1236e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 22/50
68/68 [==============================] - 0s 5ms/step - loss: 9.9636e-04 - accuracy: 1.0000 - val_loss: 4.5978e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 23/50
68/68 [==============================] - 0s 5ms/step - loss: 6.7007e-04 - accuracy: 1.0000 - val_loss: 3.9095e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 24/50
68/68 [==============================] - 0s 5ms/step - loss: 3.5613e-04 - accuracy: 1.0000 - val_loss: 3.9610e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 25/50
68/68 [==============================] - 0s 6ms/step - loss: 4.2811e-04 - accuracy: 1.0000 - val_loss: 3.3714e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 26/50
68/68 [==============================] - 0s 6ms/step - loss: 3.0813e-04 - accuracy: 1.0000 - val_loss: 3.7528e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 27/50
68/68 [==============================] - 0s 5ms/step - loss: 5.1544e-04 - accuracy: 1.0000 - val_loss: 4.0429e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 28/50
68/68 [==============================] - 0s 7ms/step - loss: 3.9930e-04 - accuracy: 1.0000 - val_loss: 3.3882e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 29/50
68/68 [==============================] - 0s 6ms/step - loss: 1.7482e-04 - accuracy: 1.0000 - val_loss: 2.2123e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 30/50
68/68 [==============================] - 0s 5ms/step - loss: 2.3158e-04 - accuracy: 1.0000 - val_loss: 2.2260e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 31/50
68/68 [==============================] - 0s 7ms/step - loss: 1.9564e-04 - accuracy: 1.0000 - val_loss: 2.0184e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 32/50
68/68 [==============================] - 0s 6ms/step - loss: 2.4279e-04 - accuracy: 1.0000 - val_loss: 1.9479e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 33/50
68/68 [==============================] - 0s 5ms/step - loss: 1.7506e-04 - accuracy: 1.0000 - val_loss: 2.1492e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 34/50
68/68 [==============================] - 0s 4ms/step - loss: 6.2078e-04 - accuracy: 0.9997 - val_loss: 1.7919e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 35/50
68/68 [==============================] - 0s 6ms/step - loss: 2.6563e-04 - accuracy: 1.0000 - val_loss: 1.5277e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 36/50
68/68 [==============================] - 0s 4ms/step - loss: 8.7712e-04 - accuracy: 1.0000 - val_loss: 1.9444e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 37/50
68/68 [==============================] - 0s 5ms/step - loss: 2.5617e-04 - accuracy: 1.0000 - val_loss: 2.9077e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 38/50
68/68 [==============================] - 0s 4ms/step - loss: 2.3949e-04 - accuracy: 1.0000 - val_loss: 1.7383e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 39/50
68/68 [==============================] - 0s 4ms/step - loss: 1.3834e-04 - accuracy: 1.0000 - val_loss: 1.5127e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 40/50
68/68 [==============================] - 0s 5ms/step - loss: 1.7125e-04 - accuracy: 1.0000 - val_loss: 2.1180e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 41/50
68/68 [==============================] - 0s 5ms/step - loss: 2.4227e-04 - accuracy: 1.0000 - val_loss: 1.3420e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 42/50
68/68 [==============================] - 0s 4ms/step - loss: 1.0846e-04 - accuracy: 1.0000 - val_loss: 1.2379e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 43/50
68/68 [==============================] - 0s 4ms/step - loss: 9.1086e-05 - accuracy: 1.0000 - val_loss: 1.1907e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 44/50
68/68 [==============================] - 0s 5ms/step - loss: 1.0132e-04 - accuracy: 1.0000 - val_loss: 1.0574e-04 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 45/50
68/68 [==============================] - 0s 4ms/step - loss: 8.9067e-05 - accuracy: 1.0000 - val_loss: 9.8274e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 46/50
68/68 [==============================] - 0s 5ms/step - loss: 7.7990e-05 - accuracy: 1.0000 - val_loss: 7.8939e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 47/50
68/68 [==============================] - 0s 5ms/step - loss: 7.0133e-04 - accuracy: 1.0000 - val_loss: 6.7715e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 48/50
68/68 [==============================] - 0s 5ms/step - loss: 1.3046e-04 - accuracy: 1.0000 - val_loss: 6.9910e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 49/50
68/68 [==============================] - 0s 5ms/step - loss: 8.6249e-05 - accuracy: 1.0000 - val_loss: 7.4691e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000
Epoch 50/50
68/68 [==============================] - 0s 6ms/step - loss: 6.3863e-05 - accuracy: 1.0000 - val_loss: 7.6789e-05 - val_accuracy: 1.0000
WARNING:tensorflow:Warning: could not reach RemoteMonitor root server at http://localhost:9000

Pipeline(steps=[('imputer', SimpleImputer()), ('scaler', StandardScaler()),
                ('model',
                 <tensorflow.python.keras.wrappers.scikit_learn.KerasClassifier object at 0x7ff9b744b1c0>)])

# Log CSV

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
plt.suptitle('Proceso de entrenamiento')
history = pd.read_csv(log_csv)
history[['accuracy', 'val_accuracy']].plot(ylabel='Accuracy', ax=ax1, xlabel='Epoch')
history[['loss', 'val_loss']].plot(ylabel='Loss', ax=ax2, xlabel='Epoch')
plt.tight_layout()
plt.show()

classifier_neigh = KNeighborsClassifier(n_neighbors=3)
classifier_neigh.fit(X=X_train, y=y_train)

KNeighborsClassifier(n_neighbors=3)

Predicciones

Predicciones en la imagen

result_keras = stack[bandsout].predict(estimator=pipeline_keras, dtype='uint8', nodata=255)
result_keras_probs = stack[bandsout].predict_proba(estimator=pipeline_keras)

result_neigh = stack[bandsout].predict(estimator=classifier_neigh, dtype='uint8', nodata=255)
result_neigh_probs = stack[bandsout].predict_proba(estimator=classifier_neigh)

588/588 [==============================] - 2s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 7))

ep.plot_rgb(rgb843, extent=extent, ax=ax1, stretch=True, title="Imagen original")
manzana.boundary.plot(ax=ax1, color='white', alpha=0.5,)

ep.plot_bands(result_keras.read(), extent=extent, cmap="RdYlGn", cbar=False, 
              title='Resultado de clasificación', ax=ax2)
manzana.boundary.plot(ax=ax2, color='white', alpha=0.5)
plt.tight_layout()
plt.show()

# Almacenamiento de los resultados de la imagen 

result_keras.write(imgkeras_classification, nodata=255)
result_keras_probs.write(imgkeras_probability, nodata=255)

result_neigh.write(imgkneig_classification, nodata=255)
result_neigh_probs.write(imgkneig_probability, nodata=255)

<pyspatialml.raster.Raster at 0x7ff9ac480b50>

Evaluaciones

Métricas de evaluación de la clasificación

# Calcular la evaluación de exactitud temática

def calculate_binary_class_scores(y_true, y_pred):
   accuracy = accuracy_score(y_true, y_pred.astype('int64'))
   precision = precision_score(y_true, y_pred.astype('int64'))
   recall = recall_score(y_true, y_pred.astype('int64'))
   return accuracy, precision, recall

# Evaluación Test clasificación de Keras

y_pred = pipeline_keras.predict(X_test)
test_accuracy_score2, test_precision_score2, test_recall_score2 = calculate_binary_class_scores(y_test, y_pred)

print('Keras Test Data Accuracy (%) = ', round(test_accuracy_score2*100,2))
print('Keras Test Data Precision (%) = ', round(test_precision_score2*100,2))
print('Keras Test Data Recall (%) = ', round(test_recall_score2*100,2))

cm = confusion_matrix(y_test, y_pred)
print("Keras accuracy is {}%".format(((cm[0][0] + cm[1][1])/y_test.shape[0])*100))
sns.heatmap(cm,annot=True, fmt='g')

22/22 [==============================] - 0s 3ms/step
Keras Test Data Accuracy (%) =  100.0
Keras Test Data Precision (%) =  100.0
Keras Test Data Recall (%) =  100.0
Keras accuracy is 100.0%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9c61411c0>

# Evaluación Test clasificación de KNeighbors

y_pred = classifier_neigh.predict(X_test)
test_accuracy_score2, test_precision_score2, test_recall_score2 = calculate_binary_class_scores(y_test, y_pred)

print('KNeighbors Test Data Accuracy (%) = ', round(test_accuracy_score2*100,2))
print('KNeighbors Test Data Precision (%) = ', round(test_precision_score2*100,2))
print('KNeighbors Test Data Recall (%) = ', round(test_recall_score2*100,2))

cm = confusion_matrix(y_test, y_pred)
print("KNeighbors Test accuracy is {}%".format(((cm[0][0] + cm[1][1])/y_test.shape[0])*100))
sns.heatmap(cm,annot=True, fmt='g')

KNeighbors Test Data Accuracy (%) =  100.0
KNeighbors Test Data Precision (%) =  100.0
KNeighbors Test Data Recall (%) =  100.0
KNeighbors Test accuracy is 100.0%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9b771e2b0>

# Evaluación Train clasificación de calidad Keras

y_pred = pipeline_keras.predict(X_train)
train_accuracy_score2, train_precision_score2, train_recall_score2 = calculate_binary_class_scores(y_train, y_pred)

print('Keras Train Data Accuracy (%) = ', round(train_accuracy_score2*100,2))
print('Keras Train Data Precision (%) = ', round(train_precision_score2*100,2))
print('Keras Train Data Recall (%) = ', round(train_recall_score2*100,2))

cm = confusion_matrix(y_train, y_pred)
print("Keras Train accuracy is {}%".format(((cm[0][0] + cm[1][1])/y_train.shape[0])*100))
sns.heatmap(cm,annot=True, fmt='g')

85/85 [==============================] - 0s 3ms/step
Keras Train Data Accuracy (%) =  100.0
Keras Train Data Precision (%) =  100.0
Keras Train Data Recall (%) =  100.0
Keras Train accuracy is 100.0%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9ac3d8c70>

# Evaluación Train clasificación de calidad KNeighbirs

y_pred = classifier_neigh.predict(X_train)
train_accuracy_score2, train_precision_score2, train_recall_score2 = calculate_binary_class_scores(y_train, y_pred)

print('KNeighbors Train Data Accuracy (%) = ', round(train_accuracy_score2*100,2))
print('KNeighbors Train Data Precision (%) = ', round(train_precision_score2*100,2))
print('KNeighbors Train Data Recall (%) = ', round(train_recall_score2*100,2))

cm = confusion_matrix(y_train, y_pred)
print("KNeighbors Train accuracy is {}%".format(((cm[0][0] + cm[1][1])/y_train.shape[0])*100))
sns.heatmap(cm,annot=True, fmt='g')

KNeighbors Train Data Accuracy (%) =  99.96
KNeighbors Train Data Precision (%) =  99.9
KNeighbors Train Data Recall (%) =  100.0
KNeighbors Train accuracy is 99.96290801186943%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9ac469bb0>

Implmenetación

Implementación dle modelo

Para imágenes SENTINEL se recomiendan los siguientes índices de vegetación:

Normalized Difference Vegetation Index (NDVI)

$$NDVI = (B8 – B4) / (B8 + B4)$$

Green Normalized Difference Vegetation Index (GNDVI)

$$GNDVI = (B8 – B3) / (B8 + B3)$$

Enhanced Vegetation Index (EVI)

$$EVI = 2.5 * (B8 – B4) / ((B8 + 6 * B4 – 7.5 * B2) + 1))$$

Advanced Vegetation Index (AVI)

$$AVI = [B8 * (1 – B4)*(B8 – B4)]^1/3$$

ndvi_file = os.path.join(path_results, f"ndvi_{img_name}-{timename}.tif")
evi_file = os.path.join(path_results, f"evi_{img_name}-{timename}.tif")

ndvi = (stack.B8.read().astype(float) - stack.B4.read().astype(float))/(
    stack.B8.read().astype(float) + stack.B4.read().astype(float)
)

## Dado que los datos de reflectancia vienen reescalados por un factor de 10.000

evi = 2.5 * (stack.B8.read()*0.0001- stack.B4.read()*0.0001)/(
    (stack.B8.read()*0.0001 + 6 * (stack.B4.read()*0.0001) 
    -  7.5 * (stack.B4.read()*0.0001)) + 1
)

raster_ndvi = Raster(np.where(ndvi >= 0.3, 1, 0).astype('uint8'), 
                     transform=result_keras.transform, crs=result_keras.crs,
                     mode='w+', in_memory=False)

raster_evi = Raster(np.where(evi >= 0.3, 1, 0).astype('uint8'), 
                    transform=result_keras.transform, crs=result_keras.crs, 
                    mode='w+', in_memory=False)

raster_ndvi.write(ndvi_file)
raster_evi.write(evi_file)

<pyspatialml.raster.Raster at 0x7ff9ac2a79a0>

# Despliege de los resultados de clasificación y NDVI

fig, axes = plt.subplots(3, 2, figsize=(8, 10))
# Plot verdadero color
ep.plot_rgb(rgb432, extent=extent, ax=axes[0,0], stretch=True, title='Verdadero color')
manzana.boundary.plot(ax=axes[0,0], color='white', alpha=0.5, linewidth=1)
# Plot falso color
ep.plot_rgb(rgb843, extent=extent, ax=axes[0,1], stretch=True, title='Falso color')
manzana.boundary.plot(ax=axes[0,1], color='white', alpha=0.5, linewidth=1)
# Plot Class Keras
ep.plot_bands(result_keras.read(), extent=extent, ax=axes[1,0], cmap='RdYlGn', 
              alpha=0.8, title="Clasificación ANN", cbar=False)
manzana.boundary.plot(ax=axes[1,0], color='white', alpha=0.5, linewidth=1)
# Plot class KNeigh
# Plot Class Keras
ep.plot_bands(result_neigh.read(), extent=extent, ax=axes[1,1], cmap='RdYlGn', 
              alpha=0.8, title="Clasificación KNeighbors", cbar=False)
manzana.boundary.plot(ax=axes[1,1], color='white', alpha=0.5, linewidth=1)
# Plot ndvi
ep.plot_bands(raster_ndvi.read(), extent=extent, ax=axes[2,0], cmap='RdYlGn', 
              alpha=0.8, title="NDVI", cbar=False)
manzana.boundary.plot(ax=axes[2,0], color='white', alpha=0.5, linewidth=1)
# Plot EVI
ep.plot_bands(raster_evi.read(), extent=extent, ax=axes[2,1], cmap='RdYlGn', 
              alpha=0.8, title="EVI", cbar=False)
manzana.boundary.plot(ax=axes[2,1], color='white', alpha=0.5, linewidth=1)
plt.tight_layout()
plt.show()

df_ndvi = raster_ndvi.extract_vector(training)
df_ndvi.rename(columns={df_ndvi.columns[0]: 'id'}, inplace=True)

df_evi = raster_evi.extract_vector(training)
df_evi.rename(columns={df_evi.columns[0]: 'id'}, inplace=True)

# Evaluación Test de NDVI 

train_accuracy_score2, train_precision_score2, train_recall_score2 = calculate_binary_class_scores(y, df_ndvi['id'].values)

print('NDVI Test Data Accuracy (%) = ', round(train_accuracy_score2*100,2))
print('NDVI Test Data Precision (%) = ', round(train_precision_score2*100,2))
print('NDVI Test Data Recall (%) = ', round(train_recall_score2*100,2))

cm = confusion_matrix(y, df_ndvi['id'].values)
print("NDVI accuracy is {}%".format(((cm[0][0] + cm[1][1])/y.shape[0])*100))
sns.heatmap(cm,annot=True, fmt='g') 

NDVI Test Data Accuracy (%) =  62.85
NDVI Test Data Precision (%) =  50.35
NDVI Test Data Recall (%) =  51.7
NDVI accuracy is 62.848664688427306%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9ac024af0>

# Evaluación Test de EVI 

train_accuracy_score2, train_precision_score2, train_recall_score2 = calculate_binary_class_scores(y, df_evi['id'].values)

print('EVI Test Data Accuracy (%) = ', round(train_accuracy_score2*100,2))
print('EVI Test Data Precision (%) = ', round(train_precision_score2*100,2))
print('EVI Test Data Recall (%) = ', round(train_recall_score2*100,2))

cm = confusion_matrix(y, df_evi['id'].values)
print("EVI accuracy is {}%".format(((cm[0][0] + cm[1][1])/y.shape[0])*100))

sns.heatmap(cm,annot=True, fmt='g') 

EVI Test Data Accuracy (%) =  62.34
EVI Test Data Precision (%) =  49.68
EVI Test Data Recall (%) =  49.41
EVI accuracy is 62.34421364985163%

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9815f3e50>

# Realizar la clasificación Keras para todas las imágenes

list_sources = list()
list_results = list()
list_titles = list()

for img in glob.glob(os.path.join(path_sources, '*.tif')):    
    imgname = os.path.basename(img)
    namepart = imgname.split('T')[0]
    imgclass = os.path.join(path_results, f"class_{namepart}-{timename}.tif")
    imgprob = os.path.join(path_results, f"prob_{namepart}-{timename}.tif")
    list_sources.append(img)
    list_results.append(imgclass)
    list_titles.append(namepart)

    newstack = Raster(img)
    [newstack.rename({name: bandsout[idx]}, in_place=True) for idx, name in enumerate(newstack.names)]
    newresult = newstack[bandsout].predict(estimator=pipeline_keras, dtype='uint8', nodata=255)
    newprob = newstack[bandsout].predict_proba(estimator=pipeline_keras)
    newresult.write(imgclass, nodata=255)
    newprob.write(imgprob, nodata=255)

588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 1ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 2s 3ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 1ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step
588/588 [==============================] - 1s 2ms/step

# Generación de todas imágenes con la clasificación  Keras

for idx in range(0, len(list_results)):  
    imgname = list_titles[idx][:8]
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 8))
    plt.axis('off')
    plt.suptitle(imgname, fontsize=25)
    
    with es.rio.open(list_sources[idx], 'r') as src:            
        ep.plot_rgb(src.read(), rgb=(6,3,2), ax=ax1,
                    stretch=True, extent=extent)
        manzana.boundary.plot(ax=ax1, color='white', alpha=0.5, linewidth=1)        
        ax1.axis('off')
        
    with es.rio.open(list_results[idx], 'r') as src:        
        ep.plot_bands(src.read(),  cmap='RdYlGn', ax=ax2, 
                      extent=extent, cbar=False)
        manzana.boundary.plot(ax=ax2, color='white', alpha=0.5, linewidth=1)
        ax2.axis('off')
            
    plt.tight_layout()
    plt.savefig(os.path.join(path_figures, f'rgb_{imgname}-{timename}.png'))
        
    plt.close()

# Generación de la animación con la secuencia imágenes con la clasificación  Keras

with imageio.get_writer(movie_classkeras, mode='I',  duration=0.8) as writer:
    for filename in sorted(glob.glob(os.path.join(path_figures, f'rgb_*-{timename}.png'))):
        print(filename)
        image = imageio.imread(filename)
        plt.tight_layout()
        writer.append_data(image)    

./tunjo/figures/rgb_20181215-20210519-094950.png
./tunjo/figures/rgb_20181225-20210519-094950.png
./tunjo/figures/rgb_20181230-20210519-094950.png
./tunjo/figures/rgb_20190213-20210519-094950.png
./tunjo/figures/rgb_20190223-20210519-094950.png
./tunjo/figures/rgb_20190325-20210519-094950.png
./tunjo/figures/rgb_20190414-20210519-094950.png
./tunjo/figures/rgb_20191021-20210519-094950.png
./tunjo/figures/rgb_20191205-20210519-094950.png
./tunjo/figures/rgb_20191210-20210519-094950.png
./tunjo/figures/rgb_20200104-20210519-094950.png
./tunjo/figures/rgb_20200109-20210519-094950.png
./tunjo/figures/rgb_20200114-20210519-094950.png
./tunjo/figures/rgb_20200213-20210519-094950.png
./tunjo/figures/rgb_20200329-20210519-094950.png
./tunjo/figures/rgb_20200413-20210519-094950.png
./tunjo/figures/rgb_20200428-20210519-094950.png
./tunjo/figures/rgb_20201129-20210519-094950.png
./tunjo/figures/rgb_20201204-20210519-094950.png
./tunjo/figures/rgb_20210113-20210519-094950.png
./tunjo/figures/rgb_20210309-20210519-094950.png

<Figure size 432x288 with 0 Axes>

reader_moviekeras = imageio.get_reader(movie_classkeras)
writer_movie = imageio.get_writer(movie_classkeras.replace('gif', 'mp4'), fps=1)

for frame in reader_moviekeras:
    writer_movie.append_data(frame)
writer_movie.close()

Análiss de cambios de coberturas vegetales multitemporal

stack_classkeras = Raster(sorted(list_results))
[stack_classkeras.rename({name: sorted(list_titles)[idx]}, in_place=True) for idx, name in enumerate(stack_classkeras.names)]

[None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None]

Se realiza la combinación de secuencias de tres imágenes (t1,t2,t3). Es decir, la alerta lo genera en las últimas tres imágenes. También se construye un histótico global

ID	Estado en el tiempo	t1 (Red)	t2 (Green)	t3 (Blue)	Color
1	Sin vegetación	0	0	0	Negro
2	Con vegetación	1	1	1	Blanco
3	Ganó vegetacion en t1 y t2	0	1	1	Cyan
4	Gano vegetación en t1	0	0	1	Azul
5	Perdió vegetación en t3	1	1	0	Amarillo
6	Perdió vegetación en t2 y t3	1	0	0	Rojo
7	Perdió vegetación en t2	0	1	0	Verde
8	Perdió en t2 y ganó en t3	1	0	1	Violeta

Para el global, se espera que los valores de celda más bajo (0) es que no ha generado ningún cambio, si cmabia a uno quiere decir que hubo un cambio de vegetación.

times_seq = list(zip(stack_classkeras.names, stack_classkeras.names[1:], 
         stack_classkeras.names[2:]))        

for idx, timeseq in enumerate(times_seq):
    if idx == 0:
        histglobal = stack_classkeras[list(timeseq)].apply(lambda x: sum(x)).read()
    else: 
        histglobal = 0.9 * histglobal + 0.1 * stack_classkeras[list(timeseq)].apply(lambda x: sum(x)).read() 
        
    ep.plot_rgb(stack_classkeras[list(timeseq)].read(), stretch=True,
                title='/'.join(timeseq))

ep.plot_bands(histglobal, cmap='inferno', title='Cambios globales')

<matplotlib.axes._subplots.AxesSubplot at 0x7ff9823cdf10>